library(ggplot2)
library(mgcv)
climate <- read.csv("https://greymonroe.github.io/PLS_206/data/climate.csv")
hemi <- read.csv("https://greymonroe.github.io/PLS_206/data/hemispheres.csv")
leaf <- read.csv("https://greymonroe.github.io/PLS_206/data/leaf_light.csv")Problem Set 4 Key: Non-linear Models and GAMs
PLS 206 — Applied Multivariate Modeling
Tip
Try the problems yourself before consulting this key. Every number below was checked in R (mgcv 1.9-4) against data/climate.csv (126 years), data/hemispheres.csv (282 rows) and data/leaf_light.csv (120 leaves) on 2026-10-18.
Part 1: Is the Warming a Straight Line?
Q1. The line rises (slope 0.0060 per year, R² = 0.683), but it is not adequate: the residuals make a U. They are mostly positive before 1900 (mean +0.14 for 1880–1900), negative from 1900 to 1980 (mean −0.10 for 1901–1920 and −0.12 for 1961–1980), and positive again after 1980 (mean +0.11). The line is too low at both ends and too high in the middle: the temperature curves upward.
m_lin <- lm(Temp ~ Year, data = climate)
summary(m_lin) # slope 0.0060, R-squared 0.683
grid <- data.frame(Year = 1880:2005)
band_lin <- cbind(grid, predict(m_lin, grid, interval = "confidence"))
ggplot() +
geom_point(data = climate, aes(Year, Temp)) +
geom_ribbon(data = band_lin, aes(Year, ymin = lwr, ymax = upr), alpha = 0.3) +
geom_line(data = band_lin, aes(Year, fit)) +
theme_classic()
climate$resid_lin <- resid(m_lin)
ggplot(climate, aes(Year, resid_lin)) + geom_point() + geom_hline(yintercept = 0) + theme_classic()Q2. Quadratic: R² = 0.801, squared term positive (7.7 × 10⁻⁵, p = 3.6 × 10⁻¹⁴). Vertex = −(−0.2918) / (2 × 7.665 × 10⁻⁵) = 1903.5, so the quadratic says the coolest year was about 1904. That is in the right neighbourhood: the coldest years in the data are 1909 (−0.48), 1904 and 1917, and the GAM in Q3 bottoms out in 1912. But the vertex is not a good estimate of the coolest period: a parabola is symmetric, so its low point is set by the steep warming after 1980 as much as by the early data. It cannot show the slight cooling from 1880 to about 1910, the bump in the 1940s, or the plateau from about 1945 to 1975.
m_quad <- lm(Temp ~ poly(Year, 2, raw = TRUE), data = climate)
summary(m_quad) # R-squared 0.801
b <- coef(m_quad)
-b[2] / (2 * b[3]) # 1903.5
band_quad <- cbind(grid, predict(m_quad, grid, interval = "confidence"))
ggplot() +
geom_point(data = climate, aes(Year, Temp)) +
geom_ribbon(data = band_quad, aes(Year, ymin = lwr, ymax = upr), alpha = 0.3) +
geom_line(data = band_quad, aes(Year, fit)) +
theme_classic()
climate[order(climate$Temp), ][1:5, c("Year", "Temp")] # 1909, 1904, 1917, 1910, 1911Part 2: A GAM for the Warming
Q3. edf = 6.80, p < 2 × 10⁻¹⁶ (F = 92.9), deviance explained = 86.2% (R² 0.683 for the line, 0.801 for the quadratic). The edf means the fitted curve uses about seven parameters’ worth of bend; an edf of 1 would be a straight line. The intercept is −0.0352, which equals the mean of Temp (−0.0352): the smooth is centred to sum to zero over the data, so the intercept carries the average level and the smooth carries only the shape.
m_gam <- gam(Temp ~ s(Year), data = climate, method = "REML")
summary(m_gam)
mean(climate$Temp) # -0.0352, the interceptQ4. At k = 10: k-index 0.66, p < 2 × 10⁻¹⁶, edf 6.8 of a maximum k′ = 9. At k = 20: edf 9.6 of 19, k-index 0.73 (p about 0.003). The curve barely changed: the two fits differ by at most 0.03 anywhere in 1880–2005, a fraction of the residual SD of 0.10. And the edf was never pressing on its ceiling (6.8 of 9, then 9.6 of 19). So k was big enough. The low k-index points to something else: the residuals are correlated in time. Neighbouring years’ residuals have a correlation of 0.35 (a warm year tends to follow a warm year), and the k-index compares neighbouring residuals, so it stays low however large k is (0.74 at k = 40). The fix is a model for correlated errors (for example gamm() with an AR(1) correlation), not a bigger k.
set.seed(206)
gam.check(m_gam) # k' = 9, edf 6.80, k-index 0.66, p < 2e-16
m_gam20 <- gam(Temp ~ s(Year, k = 20), data = climate, method = "REML")
set.seed(206)
k.check(m_gam20) # k' = 19, edf 9.59, k-index 0.73, p about 0.003
curve_ci <- function(model, label) {
p <- predict(model, grid, se.fit = TRUE)
data.frame(grid, fit = p$fit, lwr = p$fit - 1.96 * p$se.fit, upr = p$fit + 1.96 * p$se.fit, k = label)
}
bands <- rbind(curve_ci(m_gam, "k = 10"), curve_ci(m_gam20, "k = 20"))
ggplot() +
geom_point(data = climate, aes(Year, Temp), colour = "grey50") +
geom_ribbon(data = bands, aes(Year, ymin = lwr, ymax = upr, fill = k), alpha = 0.25) +
geom_line(data = bands, aes(Year, fit, colour = k)) +
theme_classic()
max(abs(predict(m_gam, grid) - predict(m_gam20, grid))) # 0.03
cor(resid(m_gam)[-1], resid(m_gam)[-nrow(climate)]) # 0.35Q5. Deviance explained: Greenhouse_gas 84.0%, Land_use 52.6%, Solar 20.4% (all three p < 0.001). Greenhouse gases track temperature best. This cannot show cause on its own because the predictors all change over the same century as the temperature: any variable with a similar trend through time will track temperature, whether or not it drives it. The models describe association with a shared time trend; showing cause needs physical understanding or experiments (here, climate physics), not a better curve.
for (v in c("Greenhouse_gas", "Land_use", "Solar")) {
m <- gam(as.formula(paste0("Temp ~ s(", v, ")")), data = climate, method = "REML")
cat(v, ": deviance explained", round(100 * summary(m)$dev.expl, 1), "%\n")
}
# 84.0, 52.6, 20.4
ggplot(climate, aes(Year, Greenhouse_gas)) + geom_line() + theme_classic()Part 3: Is the Warming Even Across the Planet?
Q6. Each smooth in a GAM is centred to sum to zero, so s(Year, by = Hemisphere) describes only each hemisphere’s shape over time, not its average level. The Hemisphere main effect carries the levels: the intercept (0.078) is the Northern mean and HemisphereSouthern (−0.051, p = 0.00015) is how much lower the Southern mean is (0.078 − 0.051 = 0.026, matching tapply()). Without it, both curves would be forced to share one average level, and the real difference between the hemispheres would have nowhere to go.
hemi$Hemisphere <- factor(hemi$Hemisphere)
ggplot(hemi, aes(Year, Temp, colour = Hemisphere)) + geom_point() + theme_classic()
m_hemi <- gam(Temp ~ Hemisphere + s(Year, by = Hemisphere), data = hemi, method = "REML")
summary(m_hemi) # HemisphereSouthern -0.051, p = 0.00015; edf 7.7 and 6.2, both p < 2e-16
tapply(hemi$Temp, hemi$Hemisphere, mean) # 0.078, 0.026Q7. Fitted anomalies (95% CI): Northern 1980 0.08 (0.03 to 0.13), 2020 1.24 (1.14 to 1.34); Southern 1980 0.20 (0.15 to 0.25), 2020 0.69 (0.60 to 0.78). The Northern Hemisphere warmed more: 1.16 vs 0.49 over the 40 years, more than twice as much. The 2020 intervals do not overlap (1.14 to 1.34 vs 0.60 to 0.78).
nd <- expand.grid(Year = c(1980, 2020), Hemisphere = levels(hemi$Hemisphere))
p <- predict(m_hemi, nd, se.fit = TRUE)
nd$fit <- p$fit
nd$lwr <- p$fit - 1.96 * p$se.fit
nd$upr <- p$fit + 1.96 * p$se.fit
nd
hgrid <- expand.grid(Year = 1880:2020, Hemisphere = levels(hemi$Hemisphere))
p <- predict(m_hemi, hgrid, se.fit = TRUE)
hgrid$fit <- p$fit
hgrid$lwr <- p$fit - 1.96 * p$se.fit
hgrid$upr <- p$fit + 1.96 * p$se.fit
ggplot() +
geom_point(data = hemi, aes(Year, Temp, colour = Hemisphere), alpha = 0.5) +
geom_ribbon(data = hgrid, aes(Year, ymin = lwr, ymax = upr, fill = Hemisphere), alpha = 0.3) +
geom_line(data = hgrid, aes(Year, fit, colour = Hemisphere)) +
theme_classic()Part 4: Light-Response Curves with nls()
Q8. Read off the plot: A at PAR = 0 averages −2.8, so Rd ≈ 2.8; the curve flattens near 21.5, so Amax ≈ 21.5 + 2.8 ≈ 24; A rises from −2.8 to 6.9 over the first 200 µmol, so alpha ≈ 9.7 / 200 ≈ 0.05. Any similar values work (the fit converged from every start we tried with Amax 15 to 30, alpha 0.02 to 0.1 and Rd 1 to 5). Estimates: Amax = 23.5, alpha = 0.060, Rd = 1.93 (all p < 0.001). The fitted Rd is smaller than the starting guess because the six dark readings are noisy; the fit uses all 60 leaves.
sun <- subset(leaf, leaf_type == "sun")
ggplot(sun, aes(PAR, A)) + geom_point() + theme_classic()
tapply(sun$A, sun$PAR, mean) # -2.84 at PAR = 0, about 21.5 at the top
m_sun <- nls(A ~ Amax * (1 - exp(-alpha * PAR / Amax)) - Rd, data = sun,
start = list(Amax = 24, alpha = 0.05, Rd = 2.8))
summary(m_sun) # Amax 23.5, alpha 0.060, Rd 1.93Q9. Amax = 23.5 (95% CI 22.5 to 24.6) µmol CO₂ m⁻² s⁻¹; alpha = 0.060 (95% CI 0.054 to 0.068) µmol CO₂ per µmol photons. Example: “Sun leaves had a maximum gross photosynthetic rate of 23.5 µmol CO₂ m⁻² s⁻¹ (95% CI 22.5 to 24.6) and an initial light-use efficiency of 0.060 µmol CO₂ per µmol photons (95% CI 0.054 to 0.068; non-linear least squares, n = 60 leaves).”
confint(m_sun)
# 2.5% 97.5%
# Amax 22.51 24.57
# alpha 0.054 0.068
# Rd 1.20 2.67Q10. Shade leaves: Amax = 9.79 (95% CI 9.00 to 10.58), alpha = 0.078 (0.062 to 0.098), Rd = 1.28 (0.58 to 1.99). Yes, they differ in Amax: the sun-leaf interval (22.5 to 24.6) is far above the shade-leaf interval (9.0 to 10.6), and sun leaves reach well over twice the maximum rate. One biological reason (any reasonable one is fine): leaves that develop in full sun build more photosynthetic capacity per unit area (thicker leaves, more of the enzymes that fix CO₂), while shade leaves are built for low light, where extra capacity would rarely be used. A plus: the alpha intervals overlap, so there is no clear difference in efficiency at low light; the two leaf types differ mainly in how much they can do in bright light.
nls_ci <- function(model, newdata, level = 0.95) {
rhs <- formula(model)[[3]]
f <- function(b) as.numeric(eval(rhs, c(as.list(newdata), as.list(b))))
b <- coef(model); fit <- f(b)
J <- sapply(seq_along(b), function(j) {
h <- 1e-6 * max(1, abs(b[j])); bj <- b; bj[j] <- bj[j] + h
(f(bj) - fit) / h
})
se <- sqrt(rowSums((J %*% vcov(model)) * J))
q <- qt(1 - (1 - level) / 2, df.residual(model))
data.frame(newdata, fit = fit, lwr = fit - q * se, upr = fit + q * se)
}
shade <- subset(leaf, leaf_type == "shade")
tapply(shade$A, shade$PAR, mean) # -1.23 at PAR = 0, about 8.5 at the top
m_shade <- nls(A ~ Amax * (1 - exp(-alpha * PAR / Amax)) - Rd, data = shade,
start = list(Amax = 9.5, alpha = 0.04, Rd = 1.2))
summary(m_shade) # Amax 9.79, alpha 0.078, Rd 1.28
confint(m_shade)
par_grid <- data.frame(PAR = seq(0, 2000, 10))
bands <- rbind(data.frame(nls_ci(m_sun, par_grid), leaf_type = "sun"),
data.frame(nls_ci(m_shade, par_grid), leaf_type = "shade"))
ggplot() +
geom_point(data = leaf, aes(PAR, A, colour = leaf_type)) +
geom_ribbon(data = bands, aes(PAR, ymin = lwr, ymax = upr, fill = leaf_type), alpha = 0.3) +
geom_line(data = bands, aes(PAR, fit, colour = leaf_type)) +
theme_classic()